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ABSTRACT 

We study the evolution of galaxy rest-frame ultraviolet (UV) colors in the epoch 4 < z < 8. We use new 
wide-field near-infrared data in the Great Observatories Origins Deep Survey - South field from the Cosmic 
Assembly Near-infrared Deep Extragalactic Legacy Survey, Hubble Ultra Deep Field 2009 and Early Release 
Science programs to select galaxies via photometric redshift measurements. Our sample consists of 2812 
candidate galaxies at z > 3.5, including 113 at z ~ 7- 8. We fit the observed spectral energy distribution (SED) 
to a suite of synthetic stellar population models, and measure the value of the UV spectral slope (/3) from 
the best-fit model spectrum. We run simulations to show that this measurement technique results in a smaller 
scatter on /3 than other methods, as well as a reduced number of galaxies with catastrophically incorrect (3 
measurements (i.e.. A/? > 1). We find that the median value of (3 evolves significantly from -1.82;!;Q g4 at z = 

4, to -2.37!jj;^^ at z = 7. Additionally, we find that faint galaxies at z = 7 have /3 = -2.68!o 24 ~2-4 after 
correcting for observational bias); this is redder than previous claims in the literature, and does not require 
"exotic" stellar populations (e.g., very-low metallicities or top-heavy initial mass functions) to explain their 
colors. This evolution can be explained by an increase in dust extinction, from low amounts at z = 7, to Ay ^ 
0.5 mag at z = 4. The timescale for this increase is consistent with low-mass AGB stars forming the bulk of the 
dust. We find no significant (< la) correlation between (3 and M^^ when measuring Mjjy at a consistent rest- 
frame wavelength of 1500 A. This is particularly true at bright magnitudes, though our results do show evidence 
for a weak correlation at faint magnitudes when galaxies in the HUDF are considered separately, hinting that 
dynamic range in sample luminosities may play a role. We do find a strong correlation between /3 and the 
stellar mass at all redshifts, in that more massive galaxies exhibit redder colors. The most massive galaxies 
in our sample have similarly red colors at each redshift, implying that dust can build up quickly in massive 
galaxies, and that feedback is likely removing dust from low-mass galaxies at z > 7. Thus the stellar-mass - 
metallicity relation, previously observed up to z 3, may extend out to z = 7 - 8. 

Keywords: early universe — galaxies: evolution — galaxies: formation — galaxies: high-redshift — ultravio- 
let: galaxies 



1. INTRODUCTION 

Recent advances in technology, including the installation 
of the near-infrared (NIR) channel of the Wide Field Cam- 
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era 3 (WFC3) on board the Hubble Space Telescope (HST), 
have revolutionized our ability to study the early universe. 
Initial WFC3 programs, including the deepest NIR image 
ever taken in the Hubble Ultra Deep Field (HUDF; ID 
11563, PI lUingworth), allowed the discovery of the first 
robust samples of galaxies at z > 7 (Bouwens et al. 2010j|; 
Oesch et all 120 lOt iMcLure et al. 
Finkelstein et all 120101; lYan et al 
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.Bunker etal.. .201c ; 
In addition to the 
exciting discovery that large samples of galaxies exist at such 
early times, these data allowed us to study the physical proper- 
ties of these galaxies, at a time only ~ 500-800 Myr removed 
from the Big Bang. 

At z 7 - 8, these WFC3 NIR data are probing the rest- 
frame ultraviolet (UV) Ught from these galaxies. While the 
Spitzer Space Telescope can observe the rest-frame optical 
light, its sensitivity is not at the level to detect the typi- 
cal, sub-L* galaxies that make up most of the published 
samples (although s tacking has resulted in some detections 
dLabbe et al.ll2010bl) , one must be cautious at these redshifts 
as high equivalent width Ha and [O lll] emission can domi- 
nate fluxes in the Spi t zer bands (e.g. iFinkels t ein et all 120091 ; 
ISchaerer & de BarrosI l20Tot IShim et al.l 1201 Ih l. However, 
the rest-frame ultraviolet alone is an important diagnostic of 
g alaxy physical proper ties. 
ICalzetti et al.l ( 119941) parameterized the UV color via the 



2 



FiNKELSTEIN ET AL. 



UV spectral slope, (3, defined as cx A^, which was de- 
signed to be measured from spectra in specified wavelength 
windows. Subsequent work has found that /3 is an excellent 
tracer of the dust extinction in galaxies, as it is well correlated 
with dust emission in the far-infrared (FIR) (defined as the 



ratio of FIR -to-UV flux), at both low redshift ^Meurer et al 
T997 l[T999l) . and z - 2 (e.g.. lSeibert e"tani2002HReddv et al 
2012). 



Studying just the rest-frame UV color at high redshift can 
thus provide powerful information on the presence of dust 
in galaxies, which in turn provides insight into the chemi- 
cal enrichment of the universe. However, at such high red- 
shifts continuum spectroscopy is beyond our grasp (though 
the upcoming James Webb Space Telescope should provide 
this capability), thus /3 is traditionally e s timated via a sin- 
gle UV color (e.g. iMeurer et alj 119971 iHathi et all 12001 
Overzier et al. 2008; Bou wens et al. l l2010bHFinkelstein et alj 
2010:.Dunlop et al . 2012). 

Using a single color. iBou wens et al] (120091) measured /3 for 
a sample of ^ 1400 galaxies from 2 < z < 6. They found 
strong evolution, in that /3 changed from ~-1.5 atz~ 2 to 
^ -2 at z 6, and found also that lower luminosity galax- 
ies appeared to be bluer at a given redshift. This implies 
that the dust content in galaxies is decreasing towards earlier 
epochs. This leads to the natural qu estion of whether galax- 
ies get even bluer at higher r edshifts. iBouwens et al.l ( l2010bh 
and iFinkelstein et all ( 120101) measured /3 for galaxies at z '--^ 
7, finding that the fa inter (mp^Mw ^, 28.3) g alaxies exhibit /? 
^ -3 in both cases. IBouwens et al.l (1201 Obi) use a measured 
small uncertainty (A/3 ^ 0.2) to investigate the possibility of 
the presenc e of extremely metal-po or (or perhaps metal-free) 
stars, while iFinkelstein et al.l (120101) argue for a larger uncer- 
tainty (A/3 ~ 0.5), and concluded that the colors of these 
galaxies are consistent with stellar populations observed lo- 
cally. Recently, McLure et al. (201 1) and Dunlop et al. (2012) 
have analyzed a sample of galaxies at z 7, and found a 
variance-weighted mean value of (3 ^ -2.0. However, we note 
that this study used more rigorous selection criteria that ex- 
cluded many of the fainter galaxies, which may explain the 
apparent redder value of the UV spectral slope. In any case, 
the differing results from previous studies highlight that more 
data are clearly needed to understand how galaxy colors are 
evolving out to the highest redshifts. 

In this paper, we use recently obtained wide field WFC3 
data in conjunction with the deep HUDF data to compile a 
much larger sample of galaxies at 4 < z < 8 to make a more 
robust determination of the evolution of /3. In §2 we describe 
the data used and catalog construction, as well as our sample 
selection techniques. In §3 we discuss our new method for 
measuring the UV spectral slope /3, as well as the simulations 
we used to test our method. In §4 we discuss the evolution 
of /3 with redshift, while in §5 we examine whether /? is cor- 
related with the UV luminosity or the stellar mass, with the 
implications discussed in §6. Our conclusions are presented 
in §7. All magnitudes are in the AB system, and throughout 
we assume a concordance cosmology with Ho = 70 km s~' 
Mpc"', n,n = 0.3 and Qa = 0.7. 



2. OBSERVATIONS 

2.1. Data 

We select our sample of galaxies from a number of 
datasets. First, we use recently obtained data from the 



Cosmic Assembly Near-Infrared Deep E xtragalactic Surve>P| 
(CANDELS; Co-PIs F aber & Ferguson: lOrogin et al.ll201lF 
iKoekemoer et alj|201 Ih . which has obtained deep WFC3/IR 
imaging of the central ~50% of GOODS-S (hereafter re- 
ferred to as GOODS-S DEEP), and shallower WFC3/IR 
imaging of the southern ^25% of GOODS-S (hereafter re- 
ferred to as GOODS-S WIDE). We also use the recent 
//577WFC3 imaging data of the HUDF main field. We 
now use all of the images of this field from the completed 
HUDF09 program, which are ~ 0.2-0.3 mag deeper in the 
F105W, F125W and F160W bands in c ontrast to the ini- 
tial papers using the H UDF09 data (e.g.. lOesch et aLllIoTol: 
IBouwens et al.l l20I0al) . Finally, we use the data from the 
WFC3 Early Release Science (ERS: I Windhorst et alJlIoTlh 
program, which imaged the northern ^25% of the Great Ob- 
serva t ories Origins Deep Su rvey (GOODS iDickinson et alJ 
I2003t iGiavalisco etani2004 Chandra Deep Fiel d - S (here- 
after, G OODS-S) with WFC3. Figure 15 of iKoekemoer et alJ 
(1201 II) highlights the relative positions of these fields. Each 
of these four fields contains imaging in the FI05W, F125W 
and F160W bands (hereafter Yios, J125 and Hieo, respec- 
tively), except for the ERS, which uses the narrower F098M 
(hereafter Y09&) in place of F105W. All of these data have 
been reprocessed in the same way; details can be found in 
IKoekemoer et alj ( 1201 lb . Throughout the paper, the term "Y- 
band" will denote either Y098 or Y105, depending on the field. 

We use archival HST optical data obtained with the Ad- 
vanced Camera for Surveys in the F435W, F606W, F775W 
and F850LP fi lters (hereafter referred to as B435, Vsoe, 1775 
and Z85o) from lBeckwith et aP (2006) for the HUDF, as well 
as the GOODS ACS v2 imaging for the DEER WIDE and 
ERS fields. All of the data has been drizzled to a pixel scale 
of 0.06". The survey areas and imaging depths probed by 
these fields are summarized in Table I . 

2.2. Object Photometry 

The size of the point-spread function (PSF) varies across 
our dataset due to the long wavelength baseline covered (~ 
0.4 - 1.8 i^Lvn). Thus, prior to photometering our imaging, we 
match the PSF of each image to the Hieo-band, which has the 
largest PSF full-width at half-maximum (FWHM). First, we 
construct empirical PSFs by stacking stars identified in each 
image, subsampling to ensure accurate alignment, and then 
sampling back up to the original pixel sc ale of 0.06". We use 
the p ECONV_TOOL routine in IDL ( IVarosi & Landsmanl 
11993b to create a convolution kernel from the PSF of each 
image and the Hieo-band PSF. This kernel is then convolved 
with each image. To ensure that the PSF of the resulting im- 
ages matches that of the Hiso-band, we measure the curves-of- 
growth of stars in each image. We iterate this procedure until 
the resulting curves-of-growth are within 1 % of the Hi6o-band 
at aperture diameters larger than 0.4" (similar to the expected 
sizes of our objects) for all bandfl 

We extract photometry for sources in these images by first 
creating a detection image for each of the four fields via an 
inverse-variance weighted sum of the J125 and H160 images. 
We then use a modified version of the Source Extractor (he re- 
after SExtractor) software (v2.8.6; lBe"rtin & Arnoutil 19961) to 

' http://candels.ucolick.org 

- Except for B435 , which was only able to be matched within < 8% of the 
Hi6o-band PSF at an aperture diameter of 0.4" (the HUDF, GOODS DEEP, 
WIDE and ERS fields had B435 curves-of-growth matched to within 3%, 4%, 
6% and 8%, respectively). 
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Table 1 

Observations Summary 



Field Area B435 Vgoe 1775 Z850 ^098 >'i05 -'125 ^160 

(arcmin^) (mag) (mag) (mag) (mag) (mag) (mag) (mag) (mag) 



HUDF 5.0 29.5 29.9 29.7 29.0 — 29.2 29.5 29.5 

GOODS-SDEEPt 85.6 28.1 28.3 27.7 27.5 — 28.0 27.7 27.5 

GOODS-SWIDE 40.1 28.1 28.3 27.7 27.5 — 27.4 27.6 27.3 

GOODS-SERS 44.7 28.1 28.3 27.7 27.5 27.6 — 28.0 27.7 

Note. — The area is measured from the image area used to detect objects, which excludes noisy regions on 
the edges. The remaining columns are 5a limiting magnitudes measured in a 0.4"-diameter apertures on non-PSF 
matched images. * The CANDELS GOODS-S DEEP field data used in this analysis represent the first four (out of 
10) epochs of J|25 and H160 data. Additionally, although it is to full depth, the CANDELS-DEEP Y105 band data 
used covers only ~ 33% of the DEEP region (the remainder will be imaged in a second epoch in December 201 1). 
We note that we selected objects over the full DEEP region; thus 67% of our CANDELS-DEEP objects do not 
yet have Y105 photometry. 



perform the photometry, which adds a buffer between the 
source and the local background cell, and removes spurious 
sources associated with the distant wings of bright objects. 
We run the modified SExtractor in two-image mode, using the 
J+H image as the detection image, and then successively us- 
ing the B435, V6O6, i775, Z850, Yo98/i{)5, Ji25 and H160 images as 
the measurement image. We used the MAP_RMS weighting 
option, and we provided accurate rms maps for each image. 
We measured photometry in small elliptical apertures, using 
the PHOT_AUTOPARAMS values of 1.2, 1.7 (similar to that 
done in Finkelstein et al. 2010), as well as in a series of cir- 
cular apertures. To test the photometric uncertainties returned 
by SExtractor, we measured the noise in the unsmoothed im- 
ages in randomly placed 0.4"-diameter apertures. We found 
that this noise was consistent with the photometric errors in 
the same size apertures returned by SExtractor, highlighting 
that the SExtractor-derived uncertainties are accurate. 

To select apertures appropriate for our science goals, we ran 
simulations in the HUDF data inserting mock galaxies into the 
images and measuring their photometry in the same manner 
as described above (these simulations are described in more 
detail in §3.2). We then compare the measured colors of these 
mock galaxies to the input color in both the small elliptical 
apertures, as well as in circular apertures with diameters of 
0.4" and 0.6", as a function of magnitude. The scatter be- 
tween the input and recovered colors for both the small ellip- 
tical and the 0.4"-diameter circular apertures are nearly iden- 
tical over the entire magnitude range, while the 0.6"-diameter 
apertures begin to have significantly worse scatter in the color 
at m > 28 (with a scatter ~ 0.05 higher than the other aper- 
tures at H160 = 29 mag). Comparing the elliptical and the 0.4"- 
diameter circular apertures, the latter yield slightly better (i.e., 
0.01 mag) results at relatively bright magnitudes (mn ^ 26- 
27), while the small elliptical apertures yield slightly better 
(i.e., 0.02 mag) results at fainter magnitudes. We then turned 
to our object catalog, and investigated which aperture yields 
better signal-to-noise measurements. We found again that at 
the faintest magnitudes, the small elliptical apertures tend to 
do better, with signal-to-noise ratios ^ 10% higher at Higo = 
29 mag, thus we adopt the fluxes measured in these apertures 
for our analysis. To correct the flux in each band to the total 
flux, we obtained the Higo-band flux in a larger Kron aperture 
with PHOT_AUTOPARAMS values of 2.5, 3.5 (the standard 
SExtractor MAG_AUTO), and then use the ratio of this flux to 
the Hi6o-band flux in the smaller elliptical aperture as an aper- 
ture correction which is applied to the fluxes in each band. 



2.3. Photometric Redshifts 

We use the EAZY software package (iBrammer et al.ll2008h 
to measure the photometric redshifts for objects in each of 
our four fields. EAZY uses non-negative linear combinations 
of templates to solve for the best-fitting photometric redshift. 
EAZY also returns the distribution, from which we can 
construct the redshift probability distribution function, P{z), 
defined as P{z) oc exp(-x^/2), normalized such that the in- 
tegral of P{z) over all redshifts is unity. We used the set of 
templates provided with EAZY b ased on the PEGASE stel- 
lar po pulation synthesis models dPioc & Rocca-Volmerangd 
119971) . We also include an additional template based on the 
rest-frame ultraviolet (UV) spectrum of the young, unred- 
dened galaxy BX418 (lErb et al.ll2010l) . as it retains character- 
istics expected in high-redshift galaxies, such as strong L yg 
emissi on. EAZY assumes the IGM prescription of Madaul 
(119951) . We note that as the luminosity functions at very high 
redshift are not yet very well-constrained, we elected to not 
include magnitude priors. 

2.4. Sample Selection 

Our goal is to study the evolution of the rest-UV color from 
z = 4 - 8. To do this, we now construct five independent galaxy 
samples at each redshift of Zsampie = 4 (3.5 < z < 4.5), Zsampie = 

5 (4.5 < Z < 5.5), Zsan,ple = 6 (5.5 < Z < 6.5), zsampie = 7 (6.5 < 

z ^ 7.5) and Zsampie = 8 (7.5 < z < 8.5). We do not base these 
samples solely on the best-fit photometric redshift; rather we 
use the information available in the full P{z) function. We do 
this by using the following selection criteria: 

1) Signal-to-noise in both J125 and Hieo > 3.5 

2) primary ^ 0.7 

3) n™„..>0.25 

4) Zsampie " 1 ^ ^phol ^ Zsampie 1 

5) xLzF<10 

where Vprimary IS the integral of P{z) in the primary redshift 
peak, and P-,,,,,,,,,, is the integral of P{z) in the redshift range of 
the sample of interest (i.e., Zsampie- 0.5 to Zsampie+ 0.5). 

The first criterion ensures that the object is significantly de- 
tected in two bands, and is thus unlikely to be spurious. The 
second criterion demands that the primary solution be a dom- 
inant one, in that only < 30% of the area under the P{z) curve 
can be in a secondary solution (see §3.4 for the effects of vary- 
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Table 2 

Galaxy Sample 



Field 


Nz=4 


N,=5 


^^=6 




^^=8 


HUDF 


177 


71 


50 


26 


14 


GOODS-S DEEP 


846 


231 


89 


19 


13 


GOODS-S WIDE 


265 


68 


23 


3 


1 


GOODS-S ERS 


647 


171 


61 


32 


5 


Total Sample 


1935 


541 


223 


80 


33 



Figure 1. Left) An example of our sample selection. This plot shows the P{z) curve for an object in our z = 5 sample. The best-fit redshift for this object is 
Zphot = 4.58, and it contains > 80% of the area under the curve in the primary redshift peak. Over 50% of the area under the P(z) curve falls in the 4.5 < z < 
5.5 (i.e., z = 5) region, while only ~25% of the area falls in the z = 4 region, thus it is placed in the z = 5 sample. Right) A comparison between the spectroscopic 
and photometric redshifts for the 157 galaxies in our sample with available spectroscopy. After excluding 12 outliers (not shown), we find excellent agreement, 
with o-:/(l+z) = 0.044. 

ing this criterion). The third criterion ensures that sufficient 
area under the P{z) curve resides in the regime of a given sam- 
ple (i.e., z = 3.5 - 4.5 for the z = 4 sample), while the fourth 
criterion ensures that the best-fit photometric redshift is near 
to the sample redshift as well. In practice, the third and fourth 
criteria individually don't exclude many sources, but the com- 
bination of the two, combined with the second criterion, en- 
sures a clean sample. The final criterion ensures that EAZY is 
providing a reasonable fit; this final criterion does not exclude 
many sources. Finally, to ensure that objects only show up 
in one sample, if an object satisfies these criteria for multiple 
samples, then it is placed in the sample where the majority of 
the area under its P{z) curve resides (i.e., if Zpiwi = 6.5, and Vj 
> Ve, then this object would be placed in the z= 1 sample). 
We highlight an example of this process in Figure 1 . We em- 
phasize that at these redshifts, while the P{z) does primarily 
key off of the Lyman break, it includes information from all 
filters, allowing the construction of an accurate P{z) function. 

After running EAZY on the catalogs for each of our four 
fields, we select galaxies in each of our five redshift bins in 
each field using the above criteria. We visually inspected each 
source, to ensure that it is real, excluding interlopers such as 
diffraction spikes or oversplit parts of bright nearby galaxies. 
We also checked for duplicates when a galaxy was in an over- 
lap region between any of our four fields; in these cases, we 
kept the object in the sample that had the deeper photome- 
try. Our final sample comprises 2812 objects, 113 of which 
are in our z = 7 and 8 samples (i.e., z > 6.5). To test the 
validity of our sample selection, we compared the photomet- 
ric redshifts of objects in our sample to known spectr oscopic 



redsh i fts in GOODS-S (e.g . iGrazian et al . 2006; Wuyts et al.l 
120081; IVanzella et al.l 120081 and references therein), finding 
excellent agreement, and very few catastrophic outliers. In 
summary, we found spectroscopic redshift matches for 157 of 
our galaxies. Of these, we found that 12 were outliers, which 
we defined as having a difference between the spectroscopic 
and photometric redshift of Az > 0.5, and that the spectro- 
scopic redshift was not consistent within 3cr of the photomet- 
ric redshift. All of these 12 had Zpiwt ~ 4, with 1 1 being lower 
redshift galaxies {Zspec < 2), and one being at higher redshift 
(Zspec = 5.1). We removed these 12 outliers from our sample. 
We found that once these outliers were removed, our pho- 



Note. — A summary of the number of candidate galaxies in 
each redshift bin from each dataset. There are a total of 2812 
candidate galaxies in our sample. 

tometric redshift accuracy was quite good, with (t-/{1+z) = 
0.044. We note that only 10 galaxies in our sample with z > 
6 have spectra, thus there is clearly a need for larger spectro- 
scopic samples at the highest redshifts. The comparison be- 
tween photometric and spectroscopic redshifts for our cleaned 
sample can be found in Figure 1, and the details of our full 
sample are summarized in Table 2. For consistency, we as- 
sume the photometric redshift in the analysis below even for 
those objects with confirmed spectroscopic redshifts. 

We note that a few of the F105W visits in the HUDF are af- 
fected by persistence from a previously imaged globular clus- 
ter The affected regions are a relatively small portion of the 
image, thus we have included all visits in our final stacked 
F105W image to obtain the maximum depth possible. To see 
if any objects in our final sample are affected, we performed 
photometry on an image incorporating the subset of ~ 75% 
of the F105W exposures unaffected by strong persistence. We 
find the difference between the fluxes of objects in our sam- 
ple between the two images is small, with a scatter of ^ 12% 
for galaxies with Yios < 29 mag. We note that much of this 
scatter is likely due to the decreased depth in the persistence- 
free image. No objects had F105W fluxes which differed by 

> 3(7, and only three objects had F105W fluxes differing by 

> 2a (two of which were near to the image edge), thus we 
conclude that the persistence is not having a significant effect 
on the F105W fluxes of the HUDF objects in our sample. 
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Figure 2. Top row: Example objects from our z = 4 sample (left) and z = 7 sample (right). The top images show 3" stamps of each object in the //5r bands. The 
plots show the best-fit model spectrum, as well as the UV spectral slopes from the SED-fitting method (gray), power-law method (blue) and single-color method 
(red). The circles show the measured fluxes of the objects in each band, while the squares show the bandpass-averaged fluxes from the best-fit model. Bottom 
row: Simulated objects at z = 4 (left) and z = 7 (right), each chosen to have a photometric redshift and Huo magnitude similar to the real object in the top panel. 
In all four cases, variations in the photometry result in a single-color-derived value of /3 which is different from the other two methods (and the true value). In 
general at z = 4, the SED-fitting and power-law methods both work fairly well, but by z = 7, the SED-fitting method appears more reliable, as the power-law 
method is reliant upon a single color. 



3. ULTRAVIOLET SPECTRAL SLOPE 

The rest-frame ultraviolet colors of galaxies can be param- 
eterized by the UV spectral slope, /3, defined in lCalzetti et alJ 
([19941) as f^ oc A^. This parameter was designed to be 
measured spectrosc opically, in specific windows defined by 
iCalzetti et al.l ( 119941) to omit spectral emission and absorption 
features. However, continuum spectroscopy is extremely dif- 
ficult at z '--^ 4, and currently impossible at z 7. Thus, /? 
is commonly estimate d at high redshift from a single rest- 
frame UV color (e.g. lMeureretal.| [T997t iHathi etal.ll2008t 
iBouwens et al.ll2010bh . where /3 is proportional to a single 
color over a redshift range of Az ^ 1. However, with our 
dataset we have multiple rest-frame U V colors over our whole 
sample except for the z = 8 bin. For example, at z ~ 4, there 
are six photometric bands redward of the Lyman break, cov- 
ering rest-frame wavelengths from ^ 1000 - 3500 A. Mak- 
ing use of all of the available photometric bands should yield 
more robust estimates of (3, as it will be less susceptible to one 
outlying data point. 

3.L Measuring /3 

We measure (3 by performing spectral energy distribution 
(SED) fitting on our sample of galaxies, finding the best- 
fitting SYnthetic_stelkr_gopu^^ using updated (201 1) mod- 
els of iBruzual & Chariot! ( I2003h . We created a suite of stel- 
lar population models assuming a Salpeter initial mass func- 
tion (IMF), varying the metallicity (0.02 < Z < 1 Zq), age 
(1 Myr < t < tn), dust extinction (0 < Ay < 2 mag, us- 



ing the attenuation law of ICalzetti et al.l ( 120001) ) and star- 
formation history (SFH oc exp"'/^, with r = 10^ 10^, 10^ 
10^ 10'°, 10", -3 X 10^ -10^ and -10'° yr). Note the neg- 
ative sign of the last three values of r; these indicate rising 
star-format ion histories, which mounting observational (e.g., 
IStarketal.| [2009: Papovi ch et all 1201 It Itee et all 1201 Ih and 
theoretical (e.g.. ,Finlator et alJl201 Ih evidence indicates is the 
d ominant mode at z > 3. T his procedure is similar to that used 
in lFinkelstein et all (120101) . thus the reader is referred there for 
more details. 

We include nebular emission lines, adding hydrogen lines 
in the Lyman, Balmer, Paschen and Bracket series, where 
the line strengths are directly proportional to the number of 
ionizing photons (which is dependent on the population age). 
Metal lines are also added, where the metal line ratios are 
a function of both the number of ionizing photons and the 
metallicit y. We used t he Lyman, Balmer and metal line ra- 
tios from IInou3 ( 201 Ih. while the Paschen an d Bracket line 
ratios come from lOsterbrock & Ferlandl (120061) . All emission 
line strengths are dependent on the ionizing escape fraction, 
which for simplicity we assume is zero (an d is likely to be < 
50% in the galaxies we observed (e.g., Yaiima et al.ll201 lb ). 
Further details on the nebular line calculations will be avail- 
able in Salmon et al. (in prep). 

We find the best-fit model via minimization, and we 
measure the value of /? directly from this best-fit spectrum 
by fitting a power la w to the spectrum usi ng the wavelength 
windows specified bv ICalzetti et al.l ( 119941) . We also measure 
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Table 3 

Properties of our Galaxy Sample 



ID 


Sample RA 


Dec 






Misoo 


log Stellar Mass 




(J2000) 


(J2000) 


(68% Conf. Range) 


(68% Conf. Range) 


(68% Conf Range) 


(68% Conf. Range) 



HUDF 339 


z4 


53.154957 


-27.764288 


3.42 (3.31 


-3.73) 


-2.29 (-2.58 - 


-2.17) 


-17.50 (-17.70- 


-17.39) 


7.73 (7.30 


- 7.95) 


DEEP 16412 


z5 


53.047127 


-27.807669 


4.48 (4.22 


-4.71) 


-2.38 (-2.41 - 


-1.86) 


-19.37 (-19.41 - 


-19.20) 


8.45 ( 8.42 


-9.38) 


WIDE 13924 


z5 


53.168690 


-27.857683 


4.53 (4.39 


- 4.66) 


-2.00 (-2.28 - 


-1.83) 


-20.55 (-20.58 - 


-20.41) 


9.14(8.95 


-9.61) 


ERS 4619 


z5 


53.158905 


-27.688990 


4.46 (4.09 


-4.85) 


-2.41 (-2.68 - 


-2.07) 


-18.73 (-18.93 - 


-18.52) 


8.45 ( 7.49 


- 8.44) 


HUDF 3233 


z6 


53.178478 


-27.788225 


5.64 (5.54 


-5.71) 


-1.95 (-2.01 - 


-1.77) 


-19.77 (-19.81 - 


-19.70) 


8.88 ( 8.57 


- 9.76) 


WIDE 14210 


z6 


53.225220 


-27.853327 


5.54 (5.44 


-5.78) 


-2.39 (-2.61 - 


-1.91) 


-20.27 (-20.39 - 


-20.18) 


8.56(8.37 


- 9.62) 


ERS 20668 


z6 


53.069393 


-27.746489 


5.56 (5.44 


-5.64) 


-2.51 (-2.54- 


-2.29) 


-20.76 (-20.77 - 


-20.67) 


9.24 ( 8.47 


-9.41) 


WIDE 8994 


z7 


53.239098 


-27.889389 


6.62 (6.34 


- 6.93) 


-2.37 (-2.68 - 


-1.50) 


-20.00 (-20.18- 


-19.85) 


9.06 ( 8.25 


- 9.55) 


HUDF 1838 


z8 


53.186256 


-27.778973 


7.56(7.34 


-7.79) 


-2.43 (-2.56 - 


-2.02) 


-19.53 (-19.66- 


-19.51) 


8.18 ( 8.09 


- 9.25) 


ERS 14661 


z8 


53.147747 


-27.723196 


7.64 (4.75 


- 8.62) 


-2.45 (-2.92 - 


-1.62) 


-19.79 (-20.03 - 


-19.47) 


8.39 ( 7.86 


- 9.78) 



Note. — A catalog of our 2812 z = 4 - 8 galaxy candidates, witli tlieir derived properties. The values in parenthesis represent the 68% 
confidence range on the derived parameters. Here we show 10 representative galaxies from our four fields and five redshift bins; the full 
table is available in machine-readable format in the electronic version. 



the absolute magnitude at 1500 A (M1500) from this spectrum 
by integrating the flux within a 100 A-wide square bandpass 
centered on 1500 A. We show an example of the results from 
this procedure in the top panels of Figure 2. We then esti- 
mate the uncertainty on /3 by running a series of Monte Carlo 
simulations, varying each object's photometry by an amount 
proportional to the photometric uncertainty, and re-deriving 
a best-fit model. In each simulation, we re-measure /3 in the 
same way as above, and then we use the distribution of /3's 
to calculate the 68% confidence range. This method has sev- 
eral advantages. First, provided the choice of stellar popula- 
tion models is appropriate and one can obtain a good fit, this 
method uses all the available colors to compute an accurate 
value of 13, rather than a single color which is much more sub- 
ject to photometric outliers. Secondly, this method obtains an 
accurate value of /3 for the specific photometric redshift. The 
single color method is designed to work over Az ^ 1, while 
the 68% confidence ranges on our photometric redshifts are 
typically Az ^ 0.3, thus fitting /3 at the best-fit photometric 
redshift should yield a more accurate result. In addition, the 
uncertainty in the redshift is accounted for in the uncertainty 
on /3, as in our Monte Carlo simulations we draw a new value 
of the redshift from the galaxy's P{z) curve in each simulation 
(setting P{z < Zphot-'^-) = 0). There is a disadvantage to this 
method as well, in that one is restricted to the range of /3 al- 
lowable in the models one chooses. In particular, for our suite 
of models, the bluest value of /3 is -3.15; if a galaxy had a 
UV spectral slope bluer than this, it would not be accurately 
recovered. However, as we show below, we do not find any 
galaxies with /3 < -3 (the bluest galaxy in our sample has a 
measured value of /3 = -2.91), thus we have not yet reached 
this limiting parameter value. Future studies probing closer to 
the epoch of first star formation will need to include models 
more appropriate to those expected stellar populations. 

We reiterate that the choice of models should have no im- 
pact on the derived values of /?. Provided the model is a 
good fit to the data, two models with differing ages or star- 
formation histories may have similar UV spectral slopes, and 
thus have similar measured values of p. To be sure this is 
the case, we eliminate galaxies with poor SED fits, defined as 
reduced > 10, from our analysis below. We find this re- 
jects very few galaxies; only 9 out of our total sample of 2812 
were excluded due to this cut. Comparing the distribution of 
/3 values from galaxies excluded by this cut to the full sample. 



we see a peak at similar values (/? ~ -2), while the excluded 
galaxies do appear to have an over-abundance of very red (/3 > 
0) colors. Had we chosen > 5 for the cut, we would not 
have excluded many more anomalously red galaxies, thus we 
move forward with our original choice of > 10. In either 
case, the conclusions below are unaffected by this choice. 

We note that since we do use all photometric bands, includ- 
ing those affected by the Lyman break and Lya emission, it 
is possible that our treatment of these two properties results 
in a slight systematic offset of the results. However, given the 
general agreement between our method, and the power-law 
method (where one fits a power-law to all photometry red- 
ward of the break; see §3.2), such an offset is likely minimal. 
We check this by comparing our best-fit j3 values to those 
obtained via fitting a power-law to the photometry. We find 
these values are consistent within their scatter, with a typical 
difference between the two methods of A/3 = 0.09 ± 0.20 at 
z = 5 and A/3 = 0.25 ± 0.40 at z = 7. We have also checked 
via our simulations (see §3.2) whether significant differences 
arise when excluding the Lyman-break bands from the fit at 
z = 4 and 5; we find essentially no difference for galaxies 
with /3 > -2, while bluer galaxies tend to be measured with 
P only 0.05 bluer when excluding the band around the break. 
We conclude that our inclusion of the Lyman break-affected 
bands are not significantly affecting our results. Finally, we 
also tested how well we recover the simulated z = 4 galaxy 
values of /3 when using only the ACS photometry, which is 
analogous to the case at z = 7, where galaxies are only de- 
tected in three bands. We find that, as expected, the scatter 
in /3 increases by ~ 0.1 and 0.2 at Higo = 24 and 29, respec- 
tively, similar to that found in the z = 7 simulations below. The 
difference in the recovered versus input values of /3 (i.e., the 
bias) are less than the scatter, at A/3 < 0.05 for bright galax- 
ies, and A/3 < 0.15 for faint galaxies. We ran a similar test 
using only the B435V606i775 photometry, which is analogous to 
the case at z = 8, and found an increase in the scatter on /3 of 
0.5 at all magnitudes. The bias was larger as well, with A/3 > 
0.5 for galaxies fainter than 28th magnitude. This shows that 
the current results at z = 8 we discuss are preliminary, due to 
the single detected color available, which may be contami- 
nated by both the Lyman break and Lya emission. In Table 3 
we list the celestial coordinates and derived values of redshift, 
stellar mass, M1500 and /3 for all objects in our sample. 
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Figure 3. Results from our simulations at = 4 and j = 7 (top and bottom panels, respectively) in the HUDF. In these simulations, we create artificial galaxies 
and insert them into the data. The galaxies are then photometered, and samples at different redshifts are then chosen in the same way as done on the actual data. 
All panels show the difference between the input value of the UV spectral slope /3 and that recovered from the photoinetry as a function of input H[5o magnitude. 
The left panel of both rows shows the results when /3 is recovered during the SED fitting process, which uses all available photometri c bands. Th e middle panel 
shows the results when using only a single rest-UV color The right panel shows the results when fitting a power law as in iBouwens et alj j2012l) . The shading 
denotes the density of siinulated galaxies, while the red squares denote the mean difference between simulated and observed values of /3 in bins of Am = 1. We 
have run similar simulations in the CANDELS fields which also show that the SED fitting method minimizes the scatter in /3. 



3.2. Testing the Method 

To test our SED-fitting method of measuring /3, we have 
run a series of simulations using the data in the HUDF. In 
these simulati ons, we created m ock galaxies with the GAL- 
FIT software (iPeng et al.l 120021) . using a range of half-light 
radii, axial ratios and Sersic indices. We assumed a distribu- 
tion of Hi6o magnitudes that is Gaussian at bright magnitudes 
and log-normal at faint magnitudes (which reproduces the ob- 
served number counts, extending the trend to fainter magni- 
tudes), and we created stellar population models with a range 
of ages, metallicities and dust extinctions to derive the fluxes 
in the remaining photometric bands. The input distribution 
of half-light radii and stellar population properties were cho- 
sen such that the output radius and J-H color distributions 
were closely matched to those of the galaxy samples. These 
"mock" galaxies were then inserted into the appropriate im- 
age. As the HUDF is only a single WFC3 field, to avoid over- 
crowding in the simulations we performed twenty iterations, 
each time inserting 1000 galaxies (forcing the positions of the 
simulated galaxies to avoid those of real objects in the image 
using the SExtractor segmentation map). We performed these 
simulations separately for each redshift bin, assuming a flat 
distribution of input redshifts. 

These images were then photometered using SExtractor in 
the same manner as was done on the original data, and pho- 
tometric redshifts were then computed using LAZY. Using 
the recovered photometry and photometric redshifts of these 
mock galaxies, we selected out redshift samples in each of 
the simulations using the same five criteria outlined in §2.4. 
In this way, the recovered galaxies from the simulations best 
match the properties of the galaxies in our sample. We then 



measured the value of /3 in three ways; first, we performed 
SED fitting on these mock galaxies using the recovered pho- 
tometry in the same way as on our object samples. Sec- 
ondly, we measured / ? using a single color, using the rela- 
tions from'Hathi et al.J ([ 2008h for the z = 4 sample, that those 
fromjDunlop et al. (2 ()12l) for the z = 5 - 8 samples (using the 
Y105-J125 color forz = 5 and 6, and the J125-H160 color for z = 
7 and 8). Thirdly, we measured /3 via a power-law fit to the 
available bandpass es, using the same ban dpasses in each red- 
shift bin as done in iBouwens et al.l (120121 see also Castellano 
et al. 2012). Examples of simulated galaxies are shown in the 
bottom panels of Figure 2. 

The results from these simulations in the HUDF are shown 
in Figure 3, where the shading denotes the density of simu- 
lated galaxies. In the top row, we show the results from our z = 
4 simulations, plotting the difference between the input value 
of /3 and the recovered value from the SED-fitting method in 
the left panel, the difference in (3 from the single-color method 
in the middle panel, and that from the power-law method in 
the right-hand panel. The bottom row shows the same three 
methods for our z = 7 simulations. In each panel, we com- 
pute the median difference (i.e., bias) as well as the scatter 
in bins of Am// = 1 mag, shown as the red squares and er- 
ror bars in both panels. Two items are immediately obvious; 
1) Both the SED-fitting and power-law methods result in a 
smaller scatter at all magnitudes than the single-color method, 
and 2) the number of galaxies "catastrophically" scattered to 
A/3 > ±1 is much less when using the SED-fitting or power- 
law methods. We note that at z = 7 the power-law and single- 
color methods suffer the same drawbacks, as the power-law 
method excludes the Y-band due to the possible contribution 
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Figure 4. Left) The recovered uncertainty in /3 from tlie SED fitting method (solid lines), the single-color method (dashed lines), and the power-law method 
(dotted lines) from our simulations in the HUDF (we expect the observed trends to shift to brighter magnitudes in the shallower datasets). At all redshifts and at 
all magnitudes, measuring /3 during the SED fitting process yields a smaller scatter than using only a single color. The SED-fitting method yields comparable 
results to the power-law method at z < 6, but it results in a much smaller scatter at ; = 7 and 8. Right) The photometric bias from our simulations. The SED 
fitting method has no bias at all but the faintest magnitudes, while the single color and power-law methods suffer a red bias at bright magnitudes, transitioning to 
a blue bias at faint magnitudes. 



by Lya emission. Our SED-fitting method still makes use of 
the Y-band photometry, allowing Lya emission in the models 
to best-match the observed Y - J125 color, then excluding the 
spectral region around Lya when fitting {3. This appears to 
work well, as the scatter in the SED-fitting method at z = 7 is 
less than the other two methods, and it does not show any bias 
due to poor fits to the Lya strength. Similarly, this method 
allows us to report the first results for /3 at z = 8, which the 
other methods cannot do (though with a large scatter). 

In Figure 4 we show the results of our simulations at all 
redshifts. The left panel shows the photometric scatter, and 
it is apparent that the single-color method suffers a much 
larger scatter at all redshifts and magnitudes than the other 
two methods. Additionally, while the SED-fitting and single- 
color methods yield similar results at z = 4 - 6, the SED-fitting 
method results in much less scatter at z = 7 - 8. The right-hand 
panel shows the photometric bias. At m/rigow < 28, the SED- 
fitting method results in essentially no bias at all redshifts, 
while both the power-law and single-color methods result in 
a bias toward redder values of /3. This red bias is likely due 
to the fixed-redshift nature of these methods. iMeurer et al.l 
( 119991) show (in their Figure 3) that a correction needs to 
be made to the value of /3 depending on the redshift of the 
galaxy; they find that for most redshifts, this correction is in 
the blueward direction, meaning that the galaxy is measured 
as redder than its true value, similar to what we find in our 
simulations. Finally, we note that at z = 7 - 8, aU three meth- 
ods contain a blue bias at faint magnitudes (at least relative to 
the z < 7 baseline), similar to the selection bias discussed in 
iDunlop et al.l (120 12h . This bias originates from the fact that 
the same 2-3 colors used to measure /3 are key to the selec- 
tion of these galaxies. For example, if the J^s-band flux of a 
galaxy fluctuates upward, it will be measured to have a bluer 
13 in our sample. However, if the Hieo-band flux experiences 
such a fluctuation, any measured /3 would be redder, thus this 
galaxy is more likely to be assigned a lower photometric red- 
shift, removing it from the sample. This bias does appear to 
be somewhat smaller in the SED fitting method at the faintest 
magnitudes (A/3 ^ -0.25) versus the other two (A/3 ~ -0.5). 
Given these simulation results, we conclude that the SED fit- 
ting method produces the most reliable estimate of /3, mini- 
mizing both the photometric scatter and bias. We thus adopt 



the value of (3 obtained from the SED fitting method in our 
analysis below. 

One possible complication to our simulations is that we 
have used the same set of models (i.e., the updated mod- 
els of Bruzual & Chariot) to create the simulated galaxies 
as we used to measure their recovered value of /3, thus one 
may worry that the SED fitting method is biased to have a 
lower scatter. To test this, we used a different set of mod- 
els for the input colors of the simulated galaxies; we used the 
older 2003 version of the Bru zual & Chariot spect ral synthe- 
sis models, with a Milky Way (ICardelli et al.ll 19891) extinction 
curve, while s t ill usin g the updated synthesis models with the 
ICalzetti et al.l (l2003) extinction curve when performing the 
SED fitting on the recovered objects. We compared the recov- 
ered uncertainty in /^ as a function of magnitude from this test 
to our baseline simulations, and found the scatter increased 
by < 5% (from H160 = 24 - 29) when using the different input 
assumptions. We thus conclude that our choice of models is 
not biasing our results. 

Finally, we note that these simulations also provide a sec- 
ondary test to our photometric redshifts which is not model- 
dependent, as the templates we use with EAZY differ from 
those used to create the simulated galaxies. We compared the 
recovered photometric redshift to the input value of the red- 
shift in each of our simulations, and we find excellent agree- 
ment, with cr./(l +z) = 0.02 - 0.05. The largest uncertainty of 
0.05 occurs in the lowest redshift bin of z = 4, likely due to 
the lack of a band completely blueward of the Lyman break 
for the lowest redshift galaxies in this bin. 

3.2. L Systematic Uncertainties 

The pure photometric uncertainties on our fluxes do not ac- 
count for any other systematic uncertainties, such as zeropoint 
offsets, aperture corrections, and other uncertainties that scale 
multiplicatively with the flux. These are typically accounted 
for by adding a small systematic error term when computing 
X^; we have thus added a 5% systematic error term when per- 
forming the SED fitting. We note that while some of these ef- 
fects may be random in nature, some may preferentially make 
fluxes in a given band brighter or fainter. For the example of 
aperture corrections, this will not affect the measured value of 
/3, as we apply a constant aperture correction to all bands for a 
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Figure 5. The rest-frame UV spectral slope /3 versus the derived absolute magnitude at rest-frame 1500 A, where the different panels display the results for the 
five redshift samples we studied. The small squares show the results for individual galaxies, where the color represents the field the galaxy was discovered in 
(Red = HUDF; Green = DEEP; Blue = WIDE; Yellow = ERS). For clarity we only show the errors on the individual points in the two highest redshift bins, but 
we do plot a characteristic error at bright (M1500 = —21) and faint (M1500 = —19) mag towards the top of each panel (shifted towards the edge for clarity). The 
gray circles show the median value of /3 measured from three subsamples split by the dashed Hues, which represent the evolving values of 0.25L* and 0.75L*. 
The uncertainties on these values are the errors on the median, derived via bootstrap Monte Carlo simulations. The Hght gray shaded region denotes the range in 
/3 spanned by typical local starbursts, while the cyan region denotes NGC1705, which is one of the bluest local galaxies, and is thought to contain little dust. 

given object (as we measure fluxes on PSF-matched images); 
thus the derived color wifl stay constant. While any zeropoint 
offsets may have a different sign and amplitude, HST has a 
remarkable ph otometric stability, thus the zeropoints are ac- 
curate to <1% (iKalirai et al.ll2009^ ■ well within the additional 
5% error term we already include. 

Nonetheless, it is interesting to examine our sample to see 
how such a systematic shift in fluxes could affect the resultant 
median value of beta. We performed this test at z = 4, as it is 
where the error on the median (3 is the smallest (see §3.3), and 
hence where any systematic uncertainty could have the largest 
impact. In this test, we systematically shifted the fluxes in a 
single band (running this test separately for both i775 and J125) 
up and down by 2%, and then re-performed SED fitting to 
measure /3. We then compared the resultant values of f3 to the 
true input value known from the simulation. When shifting 
the i775-band flux systematically for all objects up or down 
by 2%, we found that the recovered f3s differed from the true 
value by The difference was even smaller, at ^0.3% 

when making the shifts to the Ji25-band. Thus, in order to 
dominate over our measured uncertainty of at z = 4, any 
systematic error would have to be of the > 5% level, much 
greater than any expected systematic uncertainty. 

3.3. Computing Median Values of j3 

Figure 5 plots the calculated values of the UV spectral slope 
P for each of our five galaxy samples from z = 4 - 8. In these 
plots, we show the values of /3 for local star-forming galax- 
ies for reference. The light gray region denote s the range of 
UV spectral slopes spanned by the templates of iKinnev et al.l 
([199^ . which represent typical local starbursts with a range 



of dust attenuation. The cyan bar denotes NGC1705, which 
is one of the bluest local star-forming galaxies, and is thought 
to hav e sub-solar metallicity and be dust free jCalzetti et al.l 
[T994I) . 

As the uncertainty on the value of /3 in an individual galaxy 
(particularly at fainter magnitudes and higher redshifts) is 
large, we compute median values of /3 in three bins, split by 
25% and 75% of the evolving value of the cha racteristic mag- 
nitude L*, using the luminos ity functions of Bouwens et al.l 
( 120071) for z = 4, 5, and 6, and lBouwens et al.l (I2OI lb for z = 7 
and 8. These median values are shown in Figure 5 by gray 
circles. We note that our chosen value of 0.25L* for the 
faintest bin is similar to an apparent magnitude at z = 7 of 
^28.3 AB mag, which is the apparent magnitude cut used in 
iFinkelstein et al.l (120101) to split the bright from faint samples. 

To calculate the uncertainty on the median values of (3, we 
ru n bootstrap Monte Carlo simulations, similar to those done 
in IFinkelstein et al.l (120101) . which quantifies the uncertainty 
in how well we can measure the mean (and is not an esti- 
mate of the intrinsic scatter in the population). In brief, we 
first account for Poisson noise (which is most relevant in the 
highest-redshift bins) by selecting A^' = N+A x \/N galaxies in 
each bin through random sampling with replacement, where 
N is the number of galaxies in the bin, and A is a random 
number drawn from a normal distribution with a mean of zero 
and a standard deviation of unity. We then account for the 
photometric noise by taking each of the galaxies in the mod- 
ified sample, and obtaining a new estimate, /?', for /3, where 
/?' = P + B X (7p, where B is a different random number with 
the same distribution as A and erg is the photometric uncer- 
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Table 4 

Median Values of the UV Spectral Slope /3 



Sample Median /3 Median /3 Median /3 Median /3 /3-Mi5oo^ Adopted* 

All Galaxies L > 0.75L* 0.25L* < L < 0.75L* L < 0.25L* Slope M* 



z = 4 
z = 5 
z = 6 
z = 7 



1 Q9+0.00 

_i gi+0.02 
^■^^-0.06 

_2.07+<'"6 
-2 37"*^-^^ 



-2 " 



1.86 „ „, 
_l Q7+().07 

'■^'-0.04 
-2.02+*'" 

-2 ni+^-^^ 

^■OJ_A 3fl 



-1.80!'';!!^ 

-2.32+''' 
-2.68™ 
-1.88^ 



0.01 ± 0.03 
0.00 ± 0.06 
-0.10 ±0.07 
-0.20 ±0.11 
-0.03 ± 0.26 



-21.06 
-20.69 
-20.29 
-20.14 
-20.10 



Note. — * We note that although a significant slope is not observed at z = 4 and 5, if we restrict the faintest magnitude bin to 
only include galaxies from the HUDF dataset, we do recover a slope of -0.07 ± 0.01 at z = 4, and -0.09 ± 0.03 at z = 5. * The 
reference for our assumed values of the characteristic magnitude M* are Bouwens et al. (2007) for z = 4, 5 and 6, and Bouwens 
et al. (2011)forz = 7and8. 



tainty in the value of /3 for that galaxy. The median value of 
f3' is then computed from this modified sample. The simu- 
lation is repeated lO'* times. In this manner, we build up a 
distribution of median /3's in each of the bins, and the 68% 
confidence range is computed from these simulated (3 values. 
This process is done for all three bins (as well as for all galax- 
ies together) in each of the five redshift samples. We note that 
the Monte Carlo portion of these simulations are necessary to 
account for the intrinsic scatter in the rest-frame UV colors 
of galaxies, which contributes to the observed scatter in /?, in 
addition to the photometric error. The latter would be the only 
source of the observed scatter if one was was observing a pop- 
ulation of galaxies with a uniform color; this is not the case at 
lower redshift, thus we do not expect it to be the case out to 
the highest redshifts. 

We list the median values of /3 in all three luminosity bins, 
as well for all galaxies in each redshift sample, in Table 4. 
We summarize the evolution of the median /3 of all galaxies 
in each sample with redshift in Figure 6. 



3.4. Contamination 

Although the comparison with existing spectroscopy indi- 
cates contamination by lower-redshift galaxies in our sample 
is low, widespread spectroscopic coverage at the high-redshift 
end of our sample is not yet available. One way we can inves- 
tigate the possible effect of contamination on our results is to 
make our selection criteria more robust, and see how the re- 
sults, for example the median /3 at each redshift, vary. We do 
this by varying the criterion of Vprimary > 0.7. If we change 
this to a much more conservative value of > 0.9, then we find 
that the median values of f3 become ^ 0.05 - 0.15 bluer, as 
the galaxies which are removed tend to be redder. 

However, this does not mean that they are all contaminants. 
Red galaxies naturally will have a P{z) which is less robust. 
A blue galaxy with a well defined break can only have one 
redshift solution, while a red galaxy, even with a well-defined 
break, will frequently have a secondary solution assuming that 
the detected break is the Balmer, rather than Lyman, break. 
However, it is likely that our sample suffers from some level 
of contamination, and this exercise has shown that the change 
in median /3 at each redshift is small relative to the uncertain- 
ties, and that the trend we find between /3 with redshift does 
not significantly change. We do note that this more conser- 
vative cut does remove 5/8 z = 8 galaxies with /3 > 1.5; thus 
if they are truly contaminants, the median /? at z = 8 may be 
bluer than we measure. 



4. DISCUSSION 

4. 1. Colors of Faint Galaxies at 7 

We first examine the bluest data point, which is the median 
value of j3 for the faint (< 0.25L*) galaxies at z = 7, as this par- 
ticular data point has re ceived much attention in the literature. 
iBouwens et al.l (1201 Obi) first published a value of the UV spec- 
tral slope for faint z ~ 7 galaxies of /3 = -3 ± 0.2, discussing 
the possibility that these very blue galaxies could host metal- 
free star formation and/or a top-heavy initial mass function, as 
such exotic populations might be necessary to produce such 
blue values of (3. In Finkelstein et al. ( 2010), we examined 
the same dataset (the first year's release of the WFC3 data in 
the HUDF) and found a similar value of [3, but with our error 
treat ment above, we deri ved a larger uncertainty = 0.5 (see 
also IWilkins et al.ll201 Ih . We thus concluded that evidence 
for exotic stellar populations was not substantiated within the 
uncertainties (it is worth noting that one can obtain a value 
of ~ -3 even with "normal" stars, as a very young stellar 
population can create UV continua this blue, even when ac- 
coun ting for nebular cont inuum emissio iS. Also using these 
data, iDunlop et alj ( 120 12h found via simulations that there is 
a bias towards artificially blue slopes for faint galaxies, from 
which they argued the very blue colors are overestimated. Our 
simulations in §3.2 on the full-depth HUDF data do not show 
as strong a bias as that found by Dunlop et al., though we do 
find a bias of ^observed < f^intrimic by ~ 0.25 for faint galaxies 
at z = 7 and 8. 

As shown in Figures 2-4, our SED-fitting method has an 
advantage over the previously published methods which use 
only a single color. Thus, with our new method and the im- 
proved dataset, it is interesting to see how the value of the 
faint bin at z = 7 has changed. As shown in Figure 5, we mea- 
sure /? = -2.68^11 24 for faint galaxies at z = 7. This measured 
value is redder than that previously measured (-3 ± 0.5), yet 
we note that it is consistent with our previously published un- 
certainties. The updated uncertainty is smaller due to a com- 
bination of an improved measurement method and the deeper 
data in the HUDF. Although our sample is also larger due to 
the inclusion of the CANDELS and ERS fields, the faint bin 
is dominated by sources discovered in the HUDF. Account- 
ing for the selection bias discussed above (see Figure 4b), the 
value of f3 may be even redder, at /? ^ -2.4. Thus, the value 

' A stellar population with a Salpeter initial mass function (IMF), Z = 0.02 
Zq , a continuous star-formation history, zero dust and an age of 1 Myr has 
/3 = -3.15, or /3 = —2.90 when accounting for nebular continuum emission. 
These values change only slightly when assuming Solar metallicity. 
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Figure 6. Evolution of the median value of the UV spectral slope /3 with redshift for all galaxies in our sample (circles). The median M1500 for these points are: 
-19.4 fe =4), -19.7 (z =5), - 1 9.9 (; =6), -19.9 (z =7 ) and -20.0 (; =8) . The s maller symbols represent recent results from Finkelstein et al. 1 2010). McLure et al] 
1201 II) . IDunlop et alj 120121) , IBouwens et al.l 120121) . ICastellano et alj 120121) and iWilkins et al., ( .201 II) . each interpolated to represent the results at the median 
value of M1500 for our sample (with the exception of the z =7 result of Finkelstein et al. 2010, which is shown at M1500 = —18.8). We find that significant evolution 
in /3 takes place from z = 7 to 4, which is likely due to increased dust extinction. We note that the point at ; = 8 is preliminary, due to the large scatter inherent in 
measuring /3 for faint galaxies from a single color We have faded this point in the figure to caution the reader about these uncertainties. 

derived here is fully consistent with the expected UV slopes 
of "normal" stellar populations. There is no evidence for ex- 
treme or exotic stellar populations. 



4.2. Evolution of P with Redshift 

As shown in Figure 6, the median value of /? decreases sig- 
nificantly with redshift from /3 -1.8 at z = 4, to /3 ^ -2.4 
at z = 7. We note that while the measured median /3 at z = 
8 appears to be redder, the small number of z = 8 galaxies 
combined with only a single detected color currently available 
for these candidate galaxies results in large uncertainties. It 
seems unlikely that galaxies get bluer from z = 4 - 7, then red- 
der again at z = 8. However, we note that if nebular continuum 
emission begins to pl ay a dominant ro le at z = 8, this could re- 
sult in redder colors (ISchaererll2002l) . If future observations 
confirm this redder color, this could also signify an approach 
to the epoch of first galaxy formation, as an entire galaxy must 
be extremely young for nebular emission to dominate. Fo- 
cussing again on z = 7, we find that galaxies in this early epoch 
have colors very similar to NGC1705, thus they are consis- 
tent with having little-to-no dust. Our derived values of /3 are 
broadly consistent with those from previous studies, which 
have found that galaxies at z ^ 6 tend to be bluer than those 
at z ^ 3 (e.g.,[ Lehnert & Bremer 2003; Stanwav et al. 2005; 
Bouwens et al. 2006; Hathi et al. 2008; Bouwens et al. 2004 
Wilkins et al.. .201 k . Bouwens et al. ,2012; iCastellano et aP 
2OI2I) . iDunlop et aU ( 20121) recently studied the colors of ro- 
bustly detected (>8cr) galaxies at z > 5 finding that the galax- 
ies at all redshifts are consistent with (5 = -2.05 ± 0. 1 (shown 
as stars in Figure 6). This is consistent within 1 cr of our mea- 
surements, though it is in tension with our z = 7 measurement. 
Howe ver, this differ e nce is likely due to the excluded sources 
in the IDunlop et al.l ( 120121) study; if we restrict ourselves to 
the brightest bin only, we would also find (3^-2 (Figure 5 
and Table 4). 

There is a significant difference between results at z = 4, 



where ICastellano et al.l (120121) have recently published an 
analysis on the colors of z = 4 galaxies using a much more 
conservative selection, using the Veoe - H160 color in place of 
Veoe - Z850- From a sample of ~ 170 galaxies which satisfy 
their selection in the HUDF and ERS fields, they compute an 
average (3 of -2.1 ± 0.08. To investigate w hether our results 
would differ significantly, we applied the ICastellano et al.l 
(I2OI2I) color selection to our sample of ~2000 z = 4 galax- 
ies, finding that ~ 1000 galaxies satisfied the selection. The 
median val ue of /3 changed from - 1.82 to -1.96; not as blue 
as found bv ICastellano et all ( 120121) . but a significant change 
nonetheless. However, even when using this more conserva- 
tive sample selection criteria, there is still a substantial drop 
in /3 from z = 4 - 7. 

As a possible explanation for the reddening of galaxy col- 
ors with redshift, we investigate a simple model where we 
assume all galaxies in our sample formed at z = 20, have 20% 
of Solar metallicity, and are forming stars at a constant rate. 
We then find how much dust extinction needs to be added to 
the population to create colors consistent with our results at 
z = 4 and 7. We find that at z = 4, Ay = 0.5 mag of extinc- 
tion is required, while only Ay < 0. 1 mag is required by z = 
7. Thus this model is consistent with very little dust having 
been formed by z = 7. In reality, both the stellar population 
age and metallicity will be changing, but the rest-fra me UV 
slope is fairly insensitive to metallicity (c.f. Bouwen s et al.l 
l2009h . and the rest-frame UV color is dominated by the most 
recently formed stars; thus if the SFH is constant or rising, age 
is not lik ely to significantly red den the slope (see also Figure 
1 1 from IBouwens et al.l (120121) ). We note also that a vary- 
ing IMF can also affect the measured value of /3, but thus far 
the observed values are consistent with a Salpeter IMF (see 
§4.1). These derived extinction values are consistent with 
those obt ained from the relation between /3 and dust extinc- 
tion from iMeurer et al.l ( 119991) . of A^ = 0.33 mag at z = 4, 
and Ay = at z = 7. Although this relation was derived at 
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z = 0, it has been shown to work at z ^ 2 for most (i.e., not 
the youngest) galaxies (e.g. lReddv et aljr2006l) . If the galax- 
ies at z = 7 are uniformly young, then it might be that their 
infrared-exces s is lower for a give n value of /? than was as- 
sumed for t he Meur er et all (119991) relation. In this case, the 
iMeurer eTal. (1999) relation would overestimate their extinc- 
tion; thus little-to-no dust in typical z = 7 galaxies in our sam- 
ple appears likely. 

If this change in color is dominated by increased extinction 
by dust, what mechanism of dust formation is responsible? 
There are two primary sources for dust formation - super- 
novae (SNe) and asymptotic giant branch (AGB) stars. Theo- 
retical mode ls have shown that a sin gle SNe can produce > 1 
M0 of dust jTodini & Ferraral[2001l) . while observations have 
given mixed results, finding very small amounts of hot dust 
a few hundred days after the explosion, but larger amoun ts of 
cold dust many years after the explosion (iGall et al.ll201 IL and 
references therein ). Given that very dusty quasars have been 
seen at z > 6 (e.g. lWang et al.ll2008l) . it is likely that SNe are 
an efficient producer of dust in the early universe. 

However, the UV continua of typical z > 6 galaxies in 
our sample appear to be relatively unattenuated. The dif- 
ference may be due to dark matter halo mass - the dusty 
quasars at high redshift inhabit the most massive halos and 
thus likely experience accelerated evolution, while even our 
"bright" subsample of galaxies have Lf/v/ < L*, and thus in- 
habit much lower mass halos. 

In the local universe, AGB stars are the primary so urce of 
dust i njection into the interstellar medium (ISM; e.g.. lGehr3 
119891) . However, at high redshift, the universe is too young for 
low-mass stars to evolve to the AGB phase; thus AGB-created 
dust can only come from stars with main-sequence lifetimes 
less than the age of galaxy formation at a given redshift. For 
example, at z = 7, again assuming star-formation begins at z = 
20, only higher-mass stars with initial masses > 3.2 M© will 
have evolved off the main sequence (the time elapsed from 
z = 20 - 7 is ^ 600 Myr for our assumed cosmology). AGB 
stars with initial masses of 2-4 Mq are the most efficient dust 
producers, while M > 5 Mq AGB stars produce little dust (e.g. 
IGall et al.l l201lf[ 

Thus our observation of an increasing dust content from z = 
7 4 is consistent with a scenario where, at least in these sub- 
L* galaxies, SNe at z > 7 have not built up much dust. As ^ 
3-5 Mq stars enter their AGB phase at z 6-7, dust buildup 
then begins in earnest. Then, in the epoch of z = 6^, when 
~ 1-3 Mq stars enter the AGB phase, dust production ramps 
up, thus by z = 4 a typical galaxy has a significant amount of 
dust. We note that it may be that SNe at z > 7 are significant 
dust producers, and that because these galaxies are low-mass 
much of the dust may be blown out of the galaxy. We discuss 
this scenario further in §6. 1 . 

We note that dust grain growth in the ISM is an alternate, 
non-stellar, method for dust formation. In this scenario, dust 
accumulates i n the ISM on seeds formed by previous su- 
pernov ae je^.. lDraine & Salpeteil97"9l:lDwek & ScalolfT980l: 
iDraind 120091) . The timescale for this process is only a few 
X 10^ yr, thus this pr ocess can account for dust formation in 
the epoch 4 < z < 8. 'Michal owski et al.l (1201 Oh find that the 
combination of dust from SNe and AGB stars are not enough 
to explain the amount of dust detected in z > 5 quasars, and 

■* We note that lower-metallicity stars will evolve to the AGB phase a few 
hundred Myr sooner, though it is also possible that galaxy formation begins a 
few hundred Myr later, at i ~ 12 - 15 (A. Karakas, private communication). 



thus conclude that grain growth in the ISM may be a con- 
tributing factor. If z > 7 SNe are providing dust seeds, then 
grain growth could be an important dust-growth mechanism 
throughout the epoch of our observations. 

Although the timescales work for the color evolution to be 
due to AGB-formed dust, this conclusion is dependent on the 
formation redshift of these galaxies. If they do form the first 
generation of stars at z ~ 15 - 20, then they should remain 
quasi-uniformly blue until significant dust is generated. In the 
AGB-dominant scenario, this is at z < 7. Thus, better statistics 
at z = 8 can help resolve this issue. If the colors of galaxies at 
z = 8 are similar to those at z = 7, then it implies that not much 
dust is forming in that epoch, thus implying that SNe at high 
redshift are not efficient dust producers. On the other hand, 
if the colors of galaxies at z = 8 are bluer than those at z = 7, 
the gradual reddening of color from 4 < z < 8 might better 
be explained by SNe, with AGB stars playing a smaller role. 
This will be possible in the near future with upcoming data 
in the HUDF (PI Ellis) which are expected to greatly increase 
the robustness of (3 measurements at z = 8. We note that the 
current data at z = 8 (although tenuous) are more consistent 
with no evolution in /3 from z = 7 - 8 than with continued 
evolution toward bluer colors. 

5. DEPENDENCE OF /3 ON MAGNITUDE AND STELLAR MASS 

5.1. Correlation of l3 with M 1500 

As is apparent in Figure 5, we do not observe a strong rela- 
tion between /3 and the UV absolute magnitude at z = 4 and 5. 
However, at z = 6 and 7, it does appear as if fainter galaxies 
have bluer values of /?. To quantify this, we fit a first-order 
polynomial through the median data points at each redshift to 
measure the slope and its associated uncertainty. These val- 
ues are listed in Table 4. We find no significant (< la) slope 
between the galaxy colors and their absolute magnitudes at 
all redshifts. This is peculiar, as a color-magnitude relation 
has pr eviou sly been observed. Fo r example, Papovich et a|] 
(12001b and iPapovich et al.l (120041) observed a strong color- 
magnitude relation at z ~3-4. However, the color used was 
between the rest-frame UV and rest-frame optical, spanning 
the 4000 A break. This color is strongly affected by age; we 
solely use colors blueward of the 4000 A break at the red- 
shifts probed by our study, thus we are much more sensitive 
to evolution in the dust extinction than in the population age. 

iLabbe et"al] (|2007) studied the rest-frame UV colors of 
galaxies out to z ~ 3, finding a color-magnitude relation at all 
redshifts using a rest-UV single color iBouwens et al. (2009!) 
performed a simi lar analysis, finding a similar relation at z ~ 
2.5 and 4, as did IWilkins et al.l (1201 lb . who used single col- 
ors to study the dust ext inction at z = 5 - 7. On the other 
hand, iReddv et al.l ( 120081) found no significant change in the 
dust reddening over the magnitude range they probed at z 
2-3, (though their sam ple was on average b righter than ours; 
i.e., closer to L*), and iDunlop et al.l (12012b also did not de- 
tect a trend between /? and the UV magnitude, though they 
restricted their sample to galaxies which had at least one 8(t 
detection. Using a po wer-law fit to the data to measure /3, 
ICaste llano et al.lT2012b only find modest evidence for a trend 
between /3 and Mieoo; with A/3 ~ 0.2 across the magnitude 
range probed by their study (M1600 = -18.5 to -21). Recently, 
iLee et al.1 ( l201 lb found evidence for redder UV spectral slopes 
for luminous 5L*) galaxies when compared to ~ L* galax- 
ies at z ^ 4 using two colors; our somewhat limited dynamical 
range in luminosity prohibits us from observing such a trend. 
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Using a similar dataset to our own, iBouwens et al.l ( 120 12h 
recently published a study of galaxy rest-UV colors at z = 4 
- 7, computing /? by fitting a power-law to the observed pho- 
tometry (and correcting for any observational biases via simu- 
lations). From our simulations discussed in §3.2, it is unlikely 
that the different method of measuring /3 would affect the re- 
lation between /3 and the Muv, yet they find strong evidence 
for such a correlation at all redshifts while we do not. How- 
ever, the difference in the way we measure Muv may have 
an effect. As we integrate the best-fit model through a 100 
A wide bandpass centered at 1500 A, we are always probing 
rest-frame 1500 A. The g eometric mean betwee n the filters 
used in the /3 derivation by IBouwens et al.l (120 12h probes dif- 
ferent parts of the rest-frame UV depending on the redshift. 
For z = 7, this probes closer to 1500 A while at z = 4, it probes 
wavelengths redder than 2000 A, far from where our measure- 
ment probed. To see if this has an effe ct, we re-created Figur e 
5 measuring Muv in the same way as IBouwens et al.l (120121) . 
finding that defining the absolute magnitude in this way al- 
lows us to recover a >3cr significance measurement of the 
P-Muv slope at z = 4 (of -0.0 4 ± 0.01; significant, but not 
as strong as the slope found bv IBouwens et al.l (12012b ). This 
is understandable, as at z = 4, our defined absolute magnitude 
sits at one end of the wavelength baseline used for f3 (close to 
1775), which may wash out any dependance between /3 and the 
absolute magnitude. For example, if a given galaxy is brighter 
in the rest-frame UV, this corresponds to a brighter ivvs-band 
flux, which, as the one of the bluest filters, would result in a 
bluer spectral slope, washing out a correlation of redder col- 
ors with brighter magnitudes. At z = 7 this is not the case (i.e., 
Mi5()o is measured in the Ji25-band, which is in the middle of 
the three filters used to derive /3), and correspondingly we do 
see an apparent dependence between (3 and M1500, though it 
is not significant. To further investigate this, we repeated our 
analysis, defining Muv to be at 2000 A rather than at 1500 
A. This is closer to the w avelength for May at z = 4 defined 
in IBouwens et al.l (l2012h . We find that this can result in a 
marginally significant slope at z = 4, of -0.04 ± 0.02. 

Thus, the method used to define the rest-frame UV abso- 
lute magnitude can affect whether a dependance between /3 
and Muv is recovered. Additionally, although the rest-frame 
wavelengths probed in these studies are similar to ours, the 
methods of deriving the colors are different for the majority 
of them, as we use multiple rest-UV colors in our SED-fitting 
method. We tested how the measured value of (3 depends on 
these different methods. We recreated our /3-magnitude plots 
(i.e. Figure 5), deriving (3 from a single color, via SED-fitting, 
and also through a power-law fit to the photometric data. We 
made these plots using both our own photo-z selected sam- 
ple, but also using only those galaxies that satisfied standard 
LBG color selection a s well (using the selection criteria from 
IBouwens et al.l (120121) ). We found that whether one uses a 
photo-z or LBG method for the selection does not affect the 
color-magnitude dependance observed. We also found that 
we can detect a color-magnitude dependance only when us- 
ing a single color Specifically, when using a single color to 
measure /3 at z = 4, we find a slope of -0.1 1 ± 0.05, or a 2a 
significance. Thus, it would appear as if the method of mea- 
suring /3 can affect whether a color-magnitude dependence is 
observed. This may be due in part to th e bias at faint mag- 
nitudes discussed bv lDunlop et al.l (120 12h . If one only uses a 
single color to derive /3, and that color is also included in the 
selection, then a bias towards bluer colors when the signal-to- 



noise is low can be introduced. 

We also investigated whether the dynamic range in M1500 
or our choice of magnitude binning is affecting our results. 
We did this by only including galaxies from the HUDF in the 
faintest bin; this changes things very little at z > 6, but it does 
effectively move the faintest bin to ~' 1 magnitude fainter at 
z = 4 and 5. We find that this results in a significantly (7a and 
3a, respectively) detected slope of -0.07 ± 0.01 at z = 4, and 
-0.09 ± 0.03 at z = 5. The emergence of a trend when restrict- 
ing the faintest bin in this way implies that if the correlation 
between /3 and Misoo is physical, it may be flatter at bright 
magnitudes. This can be further studied in the future with the 
inclusion of a larger sample of brighter galaxies from the other 
CANDELS fi elds, and more faint galaxies from the HUDF 
paraflel fields ( lOesch et al.ll2007l;lBouwens et aLllIoTlh . 

To investigate this further, we examined the values of (3 
derived in the faint bin at z = 4 when including only HUDF 
galaxies, and also when only including CANDELSh-ERS 
galaxies, finding a difference of 0.26 ± 0.10, or 2.6a. How- 
ever, we acknowledge that the z = 4 bin is unique, in that the 
lack of U-band data may result in a somewhat higher con- 
tamination in this bin, as galaxies at the low-redshift end of 
this bin will not have the Lyman limit completely blueward of 
B435. Some evidence for this was discussed in §2.4, where we 
note that the 12 contaminant galaxies identified via the dis- 
cussed spectroscopic coverage all had photometric redshifts 
in the z = 4 bin. We note that we do not expect this contam- 
ination to be large, as the 12 removed galaxies were a small 
fraction of the total number of galaxies in our z = 4 sample 
with spectroscopic coverage. One caveat is that the major- 
ity of galaxies targetted for spectroscopy are relatively bright, 
thus we do not have an emprical constraint on the contami- 
nation of the fainter galaxies in our sample, where the largest 
discrepancy in /3 between the fields is found. 

We conclude with noting that there are a variety of system- 
atic and physical reasons why a trend may or may not be ob- 
served, with a primary factor being the definition of the UV 
absolute magnitude. While we elect to use a constant rest- 
frame wavelength to better track evolution, this may lessen 
our ability to track a correlation between /? and Mf/y. 

5.2. pStellar Mass Trend 

The fact that the detectability of a dependence between (3 
and the UV absolute magnitude can depend on the defini- 
tion of the latter implies that this may not be the most robust 
method to trace the dependence of (3, as one is comparing two 
correlated quantities. We thus investigate whether comparing 
(3 to the stellar mass results in a more robust trend. While 
the stellar mass is dependent on the absolute brightness of a 
galaxy, it is not directly dependent on the galaxy color, while 
the UV absolute magnitude is dependent when one averages 
two different filters to derive it. As we have used the SED- 
fitting method to derive the values of /?, we have also derived 
the stellar masses for all the galaxies in our sample. In Figure 
7, we again plot each of the galaxies in our sample similar to 
Figure 5, only now we plot the stellar mass on the horizon- 
tal axis. Average values of (3 were computed using similar 
bootstrap Monte Carlo simulations as before; here we split 
objects into three bins; 7 <log M/M0< 8, 8 <log MIMq< 9 
and 9 <log M/M0< 10. The median values of /3 in these bins 
are given in Table 5. 

We find a > 2(7 significance dependance between (3 and the 
stellar mass at all redshifts (except at z = 8), in that more mas- 
sive galaxies are redder. However, we caution that complete- 
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Figure 7. Similar to Figure 5, except using stellar mass on the horizontal axis. We find a significant coiTelation between /3 and the stellar mass at all redshifts. 
The solid gray lines show a Hnear fit through the median /9 points. The red curve shows the 20% completeness limit for the HUDF galaxies, derived from our 
simulations (these contours will be slightly contracted for the other fields; one can estimate where they would lie by looking at the relative distribution between 
HUDF galaxies and those in the other fields). Taking the completeness into account, it appears as if the correlation seen at z = 4 - 5 is significant, as we would 
have detected faint, red galaxies if they were there. However, this is not the case at z = 6 - 8, as the completeness contour closely traces our observed trend, 
implying incompleteness is dominating this effect in the higher-redshift bins. The dashed line in the z = 4 panel shows our estimate of the contribution of the 
intrinsic correlation of /3 with the mass-to-light ratio to the overall trend. 



Table 5 

Values of the UV Spectral Slope /3 Split by Stellar Mass 



Sample 


Median /3 
log M/Mq =7-8 


Median f) 
log M/Mq = 8 - S 


Median /3 
) log M/Mq =9-10 


# Galaxies/Bin 
Low/Middle/High Mass 


/3-Stellar Mass 
Slope 


z = 4 
z = 5 
z = 6 
z = l 
z = % 


9 99+0.13 
^■^^-0.05 

-2 40+<' 22 

-2.68^-'5 
-2.503.f 


-2 0%^ °^ 
-2 l5+0-''6 
-2 20+^-^^ 
_2.42+0-3i 
9 Tc+i];!^ 

^■•'■'-0.16 


1 00+0.02 

_j 79+0.09 

1 7o+8:?^ 
"^■''^-o.io 

_ J 75+0.23 
-1 60+^^-^^ 


52/510/927 
24/227/223 
29/104/76 
17/36/22 
3/15/12 


0.17 ±0.03 
0.30 ± 0.06 
0.40 ±0.10 
0.46 ±0.10 
0.45 ± 0.37 



ness may play a strong role here, in that low-mass blue galax- 
ies are easier to detect in our dataset than low-mass red galax- 
ies. To examine if incompleteness is dominating our trend, we 
computed completeness levels from our simulations in §3.2 
as the number of recovered galaxies divided by the number 
of input galaxies in bins of stellar mass and /?. We derived 
the 20% completeness level, and we plot this in Figure 7 as 
the red curve. We note that while these simulations were only 
performed on the HUDF images, one can estimate where sim- 
ilar contours will lie for the other regions by the distribution 
of data points for objects from the different fields. 

At z = 4, it seems our observed correlation between p and 
the stellar mass is real, as we would have detected a signif- 
icant population of red (-1.5 < /3 < -1.0), low-mass (8.5 < 
log M/Mq < 9.5) galaxies had they existed. The same is true 
at z = 5. However, our completeness begins to change at z > 6, 
in that our ability to detect low-mass red galaxies is reduced 
due to the shortened wavelength baseline which we are prob- 



ing. Inspecting the completeness contours at z = 6, 7 and 8, 
we find that the shape of the contour is similar to the shape 
of the derived /3-stellar mass slope, implying that at the high- 
est redshifts, incompleteness may be dominating the observed 
trend. However, we note that as there are relatively few red, 
low-mass galaxies at z = 4, its unlikely that a large popula- 
tion exists at z = 7. Additionally, in the z = 7 panel of Figure 
7, galaxies do not tend to pile up against the completeness 
contour, further implying that low-mass red galaxies are not 
common at very high redshift. 

One caveat here is that we are only using photometry that 
covers the rest-frame UV, while photometry covering redder 
wavelengths results in more robust estimates of the stellar 
mass (by measuring the light from lower-mass stars). One 
may thus worry that any derived correlation is due to an intrin- 
sic correlation within the models between /3 and the mass-to- 
light ratio. To investigate this, we examined our suite of stel- 
lar population models, investigating the range of mass-to-light 
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9 < log U/Uq < 10 
8 < log M/Mq^9. 
7 < log></-M'S< 8 



Figure 8. Similar to Figure 6, except here we split our sample into three bins 
by stellar mass. The red circles denote the median /3 as a function of redshift 
for galaxies with 9 < log M/Mq < 10, green circles for 8 < log M/Mq 
< 9, and blue circles for 7 < log M/Mq < 8. The green shaded region 
highlights the range in /3s from the models ofjFinlator et al. (2011) for the 
middle mass range. The dotted line shows the median value of fi from these 
inodels, while the dashed line shows the median /3 from the same inodels 
with no dust extinction. The cyan bar denotes the value of /3 for NGC1705. 
This plot indicates that the trend observed in Figure 6 is largely drived by the 
lower two mass bins. The massive galaxies have a roughly constant /3 at each 
redshift, iinplying that they are able to retain the dust formed in earlier epochs 
via SNe, while the lower-inass galaxies were not. The yellow curve (and 
the right-hand vertical axis) shows the ratio between the gas accretion rate 
and the star-formation rate using the predictions from|Papovich et al. 1 201 1, 
see §6,3), Values greater than unity imply that galaxies at that redshift are 
accreting more gas than they are converting into stars, which may hinder the 
escape of Lya photons. 

ratios at a given value of /3. Investigating -1.5 > (3 > -2.5, 
which is the observed range of /3 within our sample, we find 
that at any given value of f3, there is more than an order of 
magnitude spread in the mass-to-light ratios. This is easy to 
understand, as we are allowing a wide range of star-formation 
histories and stellar population ages, along with metallicities 
and dust attenuations, thus there is a large combination of 
model parameters which c an yield a given 0, ea ch with a dif- 
ferent mass-to-light ratio. iMcLure et alJ (1201 ll) found a sim- 
ilar result, finding their allowed models spanned a mass-to- 
light range of > 50 when looking at a given UV luminosity. 
Additionally, the mass-to-light ratio distributions for models 
with f3 ^ -2.2 and those with /3 ~ -1.8 have a >50% over- 
lap. For the non-overlapping models, the bluer tend to have 
lower mass-to-light ratios. However, the difference results in 
the redder galaxies having a mass ^3 x greater than the bluer 
galaxies; this is less than the factor of 10 or higher which we 
show in Figure 8. 

Even though there is such a large spread in M/L, we in- 
vestigate how much of an impact it may have on our derived 
relation. At z = 4, where we are the most complete, we first 
took the SED-fitting derived masses, and computed the me- 
dian M/L ratio. We then applied this median M/L to the UV 
luminosity for each galaxy to rederive a stellar mass, which is 
now independent of the individual SED-fitting results. This 
should remove any intrinsic variation of the M/L with the 
models. We then measured the resultant slope between /3 and 
the stellar mass in the same way as we did on our sample, 
and found a slope of 0.10 ± 0.07. The difference between 
our observed slope and this value is thus an estimate of the 
impact of the intrinsic M/L ratios within the models, which is 
thus a slope of 0.07. This is indicated in Figure 7 by a dashed 
line. This is ^40% of our observed trend; thus it appears as if 



our observed trend is due to both an intrinsic dependance be- 
tween (3 and the M/L ratios, as well as a dependence between 
(3 and the luminosity. In any case, these trends have a physical 
cause; a very blue galaxy will be dominated by a young pop- 
ulation of massive stars, which have low mass-to-light ratios. 

The issue of correlations between /3 and the UV M/L ratio 
could potentially be alleviated by including photometry from 
the Spitzer Space Telescope Infrared Array Camera (IRAC 
iFazio et al] 120041) when measuring the stellar mass. While 
measuring IRAC fluxes for all of the galaxies in our sam- 
ple is beyond the scope of this paper, we have measured de- 
blended IR AC fluxes for our sa mple of galaxies in the HUDF 
usmg TFIT dLaidler et al.ll2007l) . To see whether the inclusion 
of IRAC fluxes results in a significant change in the stellar 
masses, and thus our observed /3-mass correlation, we exam- 
ined the /3-mass slope for HUDF galaxies only, both with and 
without inclusion of the IRAC fluxes. We found that the de- 
rived slopes were consistent within Icr, thus we conclude that 
the lack of rest-frame optical photometry is not driving the 
observed correlations. 

For an additional test as to th e validity of our de rived 
masses, we use the simulations of iFinlator et al.l ( 1201 Ih . de- 
scribed in §6.2 below. These simulations include realis- 
tic star-formation histories which yield reasonable agreement 
with the observed rest-frame UV luminosity function and the 
growth rate of the stellar-mass density at z > 6, hence they 
represent a physically-motivated test of the robustness of our 
stellar mass estimates. The star-formation histories in these 
simulations are smoothly-rising, which is reasonable given 
the tight ness of the observed star-formation-rate - stellar mass 
relation dLabbe et al.l 120 lOal) . We have performed SED fit- 
ting using the simulated photometry of these galaxies, "ob- 
serving" the galaxies, such that their photometric data points 
are varied by an amount proportional to the photometric error 
in a given band (assuming the image depths in the HUDF). 
We then compare the best-fit stellar masses using the same 
HST ACS and WFC3 bandpasses as used on our sample to 
the "known" input value of the stellar mass. We find that we 
can recover the stellar mass very well at log M/M© > 9, while 
at log M/M© < 8 we begin to do much worse. For example, 
at z = 4, the Icr spread in recovered stellar mass at an input 
mass of log M/M© = 9 is ~ 0.25 dex, increasing to ~ 0.5 dex 
at log M/Mq = 8. At higher redshift, the lack of rest-frame 
optical is less detrimental, as the age of the Universe begins to 
become a strong constra int on the amount of mass in old stars 
dFinkelstein et al.|[2010l) . Thus, except for the lowest masses 
in our sample, our SED fits with only the rest-frame U V fluxes 
do a reasonable job of recovering the average stellar masses of 
galaxies in our sample, and any uncertainty due to the lack of 
rest-frame optical observations should not affect the median- 
(3 trends we derive. 

6. IMPLICATIONS 

6.1. Dust Formation 

As we have now shown that /3 exhibits a robust dependance 
on stellar mass, it is prudent to revisit our evaluation of the 
evolution of /3 with redshift. In Figure 8 we show again (3 
versus redshift (similar to Figure 6), only now we split our 
sample into the same three mass bins, where blue, green and 
red denote the least massive, middle-mass, and most-massive 
galaxies (split again at log M/Mq = 8 and 9). While the red- 
der (3 with increasing mass at a given redshift is easily visible, 
what is more interesting here is the evolution in (3 in galaxies 
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with a constant mass. The two lowest mass bins both show a 
trend qualitatively similar to Figure 6, in the galaxies at log 
M/M0 < 9 become increasingly bluer with increasing red- 
shift. We note again that for the lowest mass galaxies at z > 6, 
it is possible that we are missing a population of low-mass red 
galaxies, but the lack of a dominant red, low-mass population 
at z = 4 - 5 indicates that such a population is unlikely to exist 
at z > 7. 

The most massive galaxies at any given redshift show a very 
different trend, in that they appear somewhat red (/3 ^ - 1.7- 
1.9) at all redshifts. Thus our inference that perhaps SNe are 
not forming dust in the galaxies in our sample cannot be true, 
as massive galaxies at z = 7 (and possibly z = 8) already appear 
to be somewhat dusty. Since at z = 7 the universe is likely too 
young for low-mass stars to have reached their AGB phase, 
the most likely explanation for this red color is dust created 
in SNe. If dust is created in SNe in the most massive galax- 
ies, then it must be created in the lower-mass galaxies as well. 
Yet, the lower-mass galaxies still appear to be dust-free (at z = 
7 - 8 for the 8 < log M/M© < 9 sample, and at z > 5 for 
the 7 < log M/M0 < 8 sample). One plausible explanation 
is that feedback from the SNe explosions is driving the dust 
out of the low-mass galaxies, as their lower gravitational po- 
tentials make it harder for them to retain gas and dust in their 
ISM against feedback. This is testable as this same mech- 
anism could cause the metallicities of low-mass galaxies to 
evolve more strongly, thus the correlation between /3 and stel- 
lar mass might indicate the slope of the mass-metallicity rela- 
tion at these high redshifts, which will be testable in the future 
with the James Webb Space Telescope {JWST). 

What does this mean for the lowest-mass galaxies? If SNe 
were the only dust production mechanism, then one might ex- 
pect low-mass galaxies to remain somewhat dust-free. How- 
ever, we see that in the two lowest mass bins galaxies do get 
redder with time. While the lowest-mass bin is the most sus- 
ceptible to incompleteness, the 8 < log M/M0 < 9 bin shows 
a very similar trend to what we saw in Figure 6 for all the 
galaxies in the sample. These results are thus still consistent 
with the idea that AGB stars are dominating the dust creation 
in lower mass galaxies (or, at least dominating the dust that is 
retained), as it is only after z = 7 that they begin to show ev- 
idence for dust extinction. The red, massive galaxies at z = 7 
show that if SN-created dust were dominating the extinction 
in the lower mass galaxies, we would see evidence for this 
at z = 7. While dust created in AGB stars would of course 
be subject to the same feedback effects as the SN-created 
dust, this dust comes > 500 Myr after the onset of the star- 
formation episode, thus any massive stars which formed near 
the low-mass (now AGB) stars will have exploded far in the 
past, while any newly formed massive stars need not be in 
close proximity. 

6.2. Comparison with Simulations 

In Figure 8 we also show the pred i cted e volution of (3 
from the simulations of iFinlator et alj (1201 lb . These sim- 
ulations were run using a custom version of the cosmo- 
logical hydrodynamic code GADGET-2 that includes recipes 
for cooling, star formation, and momentum-conserving out- 
flows. The simulated stellar continua are derived by inte- 
grating over the discrete stellar populations in each galaxy 
using the Bruzual & Chariot (2003) population synthe- 
sis models. Nebular line emission is treated following 
ISchaerer & de BarrosI (l2()lCf). and dust is added via the fore- 
ground screen model of ICalzetti et al.. (.2000i) with the nor- 



malization E(B-V ) tied to the galaxy metallicity following a 
local calibration ( IFinlator et alj|20()6h . The resulting simu- 
lations reproduce reasonably well the slope, evolution, and 
normalizatio n of the observed galaxy mass-metallic ity rela- 
tion at z < 3 (I Finl ator & Dave"2008': !Dave et al.ll201 ll) as well 
as the observed specific star formation rat es and luminosity 
functions of Lyman -break galaxies at z > 6 dPave et alj|2006l : 
IFinlator et al.ll20ll1) . The good agreement that these models 
demonstrate with existing constraints indicates that we may 
use them to understand the possible implications of our new 
observations. 

In Figure 8 we show the values of p and stellar mass from 
the simulations; to be sure there was no bias we also mea- 
sured these parameters from the simulated galaxy photometry 
in the same way as for our observed objects. We found ex- 
cellent agreement between our SED-fitting derived values of 
/3 and stellar mass and the input values from the simulated 
galaxies. These simulations have relatively few galaxies in 
our massive bin, while our observed galaxies are more rare 
in the lowest-mass bin, so we compare to the simulations in 
the middle bin of 8 < log M/M© < 9, where both the simula- 
tions and observations have a large number of galaxies. The 
simulations are denoted by the green shaded region, which de- 
notes the 68% spread in (3 values at a given redshift. We find 
good agreement between our observations and these simula- 
tions over the entire redshift range probed. We note that the 
simulations of Dayal and Ferrara (2012) also predict a depen- 
dance of [3 on the stellar mass, albeit at a lower level (A/3 ~ 
0.2 at z = 7 from log M = 7.5 ^ 9.5). 

This agreement has a number of implications. First, the 
simulations allow us to examine the galaxies with and with- 
out dust, which in turn allows us to examine how much of the 
trend of increasing /3 with redshift is driven by an increase in 
dust. We find that the fully processed simulated galaxies (i.e., 
including dust extinction) have a /3 that rises A/3 ~ 0.5 from 
z = 8 - 4 (dotted green line in Figure 8), while the underly- 
ing stellar spectrum rises by A/3 ^ 0.2 during that same time 
(dashed green line in Figure 8). Thus, while increasing age 
and metallicity is contributing to the reddening of the color, 
we confirm that dust extinction is the dominant factor in the 
model galaxies. 

As mentioned above, these simulations assume a constant 
value of extinction to metallicity, using resu lts from local 
galax ies in the Sloan Digital Sky Survey (see IFinlator et alj 
I2OO6I) . Since the evolution of /3 at a common stellar mass 
within these simulations appear to agree with our observa- 
tions, and /3 is dominated by dust, this then implies that the 
way in which galaxies suffer extinction does not change much 
over ~' 13 Gyr of cosmic time. 

Additionally, this consistency of /3 between the simulations 
and our observed galaxies imply that our observed galaxies 
may have similar metallicities to those in the simulations at 
a given redshift and stellar mass. We computed the average 
and standard deviation of the simulated galaxy mass-weighted 
stellar metallicities at each redshift for galaxies with 8 < log 
M/M0 < 9. We find that the typical galaxy in this mass range 
likely has a stellar metallicity of Z = 0.22 ± 0.06 Zq at z = 4, 
dropping to Z ~ 0. 13 ± 0.04 Zq by z = 7. These metallicities 
are low, but are well above some of the most metal-poor local 
galaxies, adding further evidence that these galaxies do not 
host Population 111 stars. 

Finally, these models include the effects of outflows, 
which reproduces the observed mass-metallicity relation at 
lower-redshifts (e.g. .Tremonti et al.. .2004; ,Erb et al. .2006h 
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by removing more gas (and metals) in lower-mass galaxies 
jFinlator et al.ll201 UlPave et al.ll201 lb . This agrees well with 
our observed dependance of /3 with mass, which at its heart 
is a mass-dust relation, and thus implies that the mechanisms 
which drive the z < 3 mass-metallicity relations exist out to 
the highest redshifts we can observe. In fact, our derived 
/3-mass slope at z = 4 (where we are the most complete) of 
0.17 ± 0.02 is consistent with th e value of 0.15 we mea- 
sure from the iFinlator et al.l (1201 lb models, further implying 
that the slope we observe is in fact a manifestation of the 
mass-metallicity relation (see also e.g., Bouwens et al. 2009; 
Bouwens et al. 2011 for similar evidence for a luminosity- 
metallicity relation at z > 4). 

This work motivates a more detailed theoretical inquiry into 
the detailed dependence of gas fraction, UV continuum color, 
and metallicity on stellar mass as well as the evolution of these 
trends. 

6.3. Probing Reionization 

Finally, we come back to the increasingly blue colors of our 
overall sample of candidate galaxies with increasing redshift. 
The evidence points to these galaxies becoming progressively 
more deficient in dust (and possibly metals) out to z = 7. Thus 
one may expect that the incidence of strong Lya emission 
from these galaxies should be increasing as well. In fact, this 
is seen out to z = 6 jStark et al.ll201(]ll201 lb . However, recent 
work has shown that the incidence of detectable Lya emission 
in galaxies decreases strongly from z = 6 to 7, with the infer- 
ence being that the volume neutral fraction of hydr ogen in the 
intergalactic medium (IGM) rises up to ^ 50% (e.g . lOno et al.l 
I2012t iPentericci et al.ll20Tll ISchenker et al.ll2012l) . However, 
it is also possible that gas within the galaxy, rather than in 
the IGM, also attenuates Lya. The increasingly blue colors of 
galaxies in our sample could indicate that these galaxies are 
becoming more gas-rich towards higher redshift, as a higher 
gas surface dens ity begets more star-formation, resulting in 
bluer colors (e.g. iKennicutBl 1 998b . 

Using observations of the star-formatio n rates, stellar 
masses and sizes of galaxies from z = 2-8, iPapovich et al.l 
(Hoil) inferred the evol ution of the gas accret ion rate with 
redshift (see Figure 4 in IPapovich et al.l (1201 lb ). In Figure 
8, we also plot their value for the ratio of the gas accretion 
rate (Mjos) to the star-formation rate (Msfr). This value is 
predicted to rise strongly with redshift; at z ^ 3, galaxies are 
turning all of the accreted gas into stars, while by z ~ 7, galax- 
ies are accreting ^ 2.5 x more gas than they can convert into 
stars. In fact, this ratio of the gas accretion rate to the star- 
formation rate rises by ^ 40% from z = 6 to z = 7. 

Therefore, one interpretation based on these model predic- 
tions and the relative blue colors of the z ^ 7 galaxies is that 
these galaxies are accreting vast reservoirs of gas. This gas is 
likely to have an attenuating effect on the Lya emission line. 
We note that gas outflows in galaxies have been invoked as 
an escape mechanism for Lya emission, as Lya photons es- 
cape from outflowing gas that is redsh ifted with respect to the 
systemic redshift of the galaxy (e.g.. IVerhamme et al.ll2008t 
iDiikstra & Wvithel2010b . However, at z 7 if the gas outflow 
rate is related to the SFR, and the gas accretion rate exceeds 
the SFR, it stands to reason that the gas inflow rates dominate 
the gas dynamics in these galaxies. The gas inflows would 
have a net blueshift with respect to the galaxy systemic red- 
shift, which would effectively block any Lya emission from 
the galaxy. Therefore, the decrease in the incidence of Lya 
in galaxies from z ^ 6 to 7 could be explained in part by an 



increase in the gas accretion rate and a decrease in the gas 
outflow rate going to higher redshift, and thus may still be 
consistent with a mostly ionized IGM. 

Finally, we note that this scenario can also explain the ob- 
servation that the Lya incidence drops mor e in faint galax- 
ies tha n bright galaxies (lOno et al. 2012) . Krumholz & Dekell 
(12012b have shown that as the SFR depends on the metallic- 
ity, galaxies in smaller dark matter halos are less able to effi- 
ciently turn gas into stars than those in large halos. In these 
smaller galaxies the suppressed SFR results in a lower metal 
production, which keeps the SFR much lower than the gas in- 
flow rate, and thus likely suppresses Lya emission. However, 
the observed mass-metallicity relation at z = implies that 
massive galaxies are more able to retain their metals, hence 
they would advance past the stage of not being able to keep up 
with inflows sooner than less-massive galaxies, which would 
cause their Lya EWs to evolve less rapidly and their metallic- 
ities (and hence dust columns) to evolve more rapidly, just as 
observed. 

Further statistics of the Lya incidence, as well as gas 
mass measurements with ALMA and metallicity measure- 
ments with JWST will allow for more detailed tests of these 
models. 



7. CONCLUSIONS 

Using the new deep, wide-area data from the CANDELS 
Multi-cycle HST program, as well as existing data from the 
HUDF and ERS programs, we have analysed the evolution 
of the UV spectral slope /3 with redshift. We utilized photo- 
metric redshifts to select galaxy samples at z = 4, 5, 6, 7 and 
8. We measured f3 with a new technique, fitting a slope to 
the best-fit model spectrum during the SED fitting process in 
the Calzetti et al. (1994) defined wavelength windows. Us- 
ing simulations, we have shown that this method determines 
/3 with smaller scatter and bias than the previously used meth- 
ods at all redshifts and magnitudes. 

We computed the median value of /3 in each redshift bin, as 
well as in subsamples split by both UV luminosity and stellar 
mass, and reached the following conclusions: 

• The median value of /3 for faint galaxies at z = 7 is /3 
= -2.68;!;q 24 (~ -2.4 after correction for observational bias). 
This value is redder than those previously reported in the lit- 
erature (/3 ~ -3), and has an uncertainty reduced by ~ 30%. 
These combine to cast further doubt upon the existence of a 
large population of metal-free stars at high redshift, as /3 = 
-2.7 is consistent with a young (< 50 Myr), dust-free popula- 
tion of 0.2 Zq stars. 

• We measure significant evolution of the median /3 for 
galaxies at all luminosities from z = 4 - 7. Galaxies at z = 
7 are again consistent with little-to-no dust, while those at z = 
4 have rest-UV colors consistent with significant dust extinc- 
tion on the order of Ay ~ 0.5 mag (depending on the exact 
formation redshift of the galaxies). If the typical galaxy in 
our sample forms its first generation of stars at z ~ 15-20, 
then this onset of dust extinction at z < 7 is commensurate 
with the evolution of M < 3.5 M0 stars off the main sequence 
and through their AGB phase. AGB stars with initial masses 
of ^2-4 Mq are the most efficient dust producers, thus this 
trend in /3 may be indicative of the increased influence of dust 
formed in the atmospheres of AGB stars. The low level of dust 
at z = 7 might indicate, at least for the typical L < L* galaxies 
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in our sample, that dust formed in SNe plays a smaller role. 

• We fit for a linear correlation between (3 and M1500 to 
determine whether we detect a significant /3 - UV absolute 
magnitude relation. At z > 6, we detect no significant relation, 
however the uncertainties on this slope are high, so a more 
robust determination of whether such a slope exists will need 
to wait for larger high-redshift samples. However, at z = 4 
and 5, we also measure slopes consistent with zero at higher 
significance, though we note that this may be dominated by 
a lack of a trend at bright magnitudes, and a possible trend 
at fainter magnitudes. Thus, the dynamic range in Muv is 
crucial, as when restricting our faint bin to HUDF galaxies 
only, we do observe a correlation. To compare our results to 
those in the literature, we have run various tests, and we find 
that the definition of the rest-frame wavelength used for Mj/v 
can determine whether a correlation is observed or not. Our 
method of using a constant rest-frame wavelength of 1500 A 
may wash out such a correlation, as we detect a significant 
correlation at z = 4 when using 2000 A to define M^uv (though 
the slope of the correlation is much less than other studies 
have found). Additionally, the method of deriving /3 can affect 
the results, as if one uses only a single color to derive /?, one 
can observe such a p-Muv dependence at the 2cr level (though 
we note some studies do acknowledge and correct for this). 
Future work examining this correlation with a larger sample 
covering a wide range of Muv, and examining all possible 
biases, will be crucial to understand the origin of any such 
correlation. 

• Comparing /3 to the stellar mass does result in a significant 
correlation, in that more massive galaxies appear redder This 
may imply that more massive galaxies experience accelerated 
evolution, as they appear red at all redshifts in our sample, 
though selection effects make this difficult to quantify at the 
highest redshifts. It thus may be that SNe are efficient dust 
producers, but only the most massive galaxies (> 10^ M©) 
can retain this dust, while lower-mass galaxies lose the dust 
in outflows. In the low mass galaxies, it is not until lower- 
mass stars pass through the AGB phase that their dust masses 
begin to build in earnest. 

• We find excellent agreement between our observed evo- 
lution o f /? with redsh i ft and the same trend in the simula- 
tions of iFinlator et al.l (1201 lb . which replicate the observed 
mass-metallicity relation at lower redshifts. These simula- 
tions assume the low-redshift value for amount of extinction 
per metallicity in a galaxy, thus this agreement implies that 
galaxies at high redshift contain a similar yield of dust-per- 
metal as at low redshift. Additionally, the slope of the (3 - 
ma ss relation we measur e at z = 4 is consistent with that from 
the IFinlator et al.l ( 1201 ll) simulations, implying that the ob- 
served evolution in the mass-metallicity relation continues out 
to at least z = 4. 

• The blue colors of our low-mass galaxies indicate that 
they should have a strong Lya emission line, yet recent results 
have shown that galaxies at z = 7 have weaker than expected 
Lya emission. Recent work indicates that galaxies are likely 
building vast reservoirs of gas with increasing redshift, which 
is plausible given the blue colors of our galaxies. This may in 
turn reduce the visibility of the Lya line without the need to 
invoke a large neutral fraction in the IGM. 

While the evolution of galaxy rest-frame UV color from z = 
4 - 7 is robust, the mechanism behind this trend is less clear 



It is likely due to dust, but the exact interplay between dust- 
formation, galaxy stellar mass, and feedback requires further 
study. More robust results at z = 8 (i.e., utilizing data from 
the upcoming treasury proposal in the HUDF; PI Ellis) may 
shed light on the issue, allowing examination of whether the 
observed trends continue into that earlier epoch. 
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